Renewable and Sustainable Energy Reviews 16 (2012) 29-36 


Contents lists available at SciVerse ScienceDirect 


Renewable and Sustainable Energy Reviews 


(ow 


ELSEVIER journal homepage: www.elsevier.com/locate/rser 


Development of a simulation model for a three-dimensional wind velocity field 
using Ténés Algeria as a case study 


Fouad Boukli Hacéne’, Miloud Tahar Abbés?, Nachida Kasbadji Merzouk‘, Larbi Loukarfi?, 
Hacene Mahmoudi?:“*, Matheus F.A. Goosen € 


a Faculty of Sciences, University Hassiba Benbouli, Bp 151, Chlef, Algeria 

> Faculty of Technology, University Hassiba Benbouli, Bp 151, Chlef, Algeria 

€ Solar Equipments Development Unit (UDES), Center for Renewable Energies Development (CDER), BP. 386, Bou-Ismail, 42415, Tipaza, Algeria 
d Alfaisal University, Riyadh, Saudi Arabia 


ARTICLE INFO ABSTRACT 

Article history: , A mathematical simulation model was developed that can determine the three-dimensional wind veloc- 
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Available online-13 September 201] edge of wind fields is crucial for predicting the dispersion of pollutants, for forecasting meteorological 


weather, for fire spread prediction and in the design and implementation of wind turbines. By means 
of a mass consistent model, an in-house program was developed to calculate the three-dimensional 
wind velocity field in the study region. The model was supported by a numerical box in which flow 
through is allowed for in the upper and lateral boundaries. The bottom boundary through which no flow 
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Wind fields simulation through occurs was determined by the topographic relief at the surface. From measured wind veloci- 
Finite difference method ties, observed values were calculated by interpolation-extrapolation. Using an optimization method, the 
Complex terrain adjusted velocities were obtained from constraints, observed velocities and the continuity equation. The 


model was verified with wind point data, the relative error did not exceed 6%. 
© 2011 Elsevier Ltd. All rights reserved. 
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1. Introduction 
The interest in wind energy production has increased consider- 
ably due to improvements in turbine technology and the growing 
need for power production from renewable sources [1-3]. One of 
* Corresponding author at: Faculty of Sciences, University Hassiba Benbouli, Bp the most important sues to deal with is the choice of an adequate 
151, Chlef, Algeria. site for a wind farm [4]. Different methodologies have been devel- 
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Nomenclature 

Xr Lagrange multiplier 

Oi variances of the observed field 

Qi Gauss precision module 

F objective function 

V Computational domain 

u,v? andw? fixed observed velocity value in x, y and z direc- 
tion 

u,vandw velocity adjusted value in x, y and z direction 

Ax, Ay, Az step in x, y and z direction 


actual measured data [5]. Mortensen and Petersen [6] for exam- 
ple have shown that numerical models are suitable for simulating 
wind fields. Such models can indicate the best sites for installation 
of wind turbines. 

Due to a boost in the world’s energy consumption as a result 
of population growth and industrialization, renewable energy (i.e. 
solar, wind, geothermal and wave) has increased in importance [7]. 
A study by Ludwig and Byrd [8] has shown that meteorological and 
air pollution problems frequently require the interpolation of com- 
plete, three-dimensional wind field flows from a limited number of 
observations at specific location. Modeling atmospheric flow is dif- 
ficult, especially when simulating wind fields over irregular terrain 
[9]. 

Different types of wind field models have been developed in 
the last few years [10]. They can be classified into two main mod- 
els: diagnostic (kinetic) and prognostic (dynamic). Kinetic models 
are those capable of reconstructing a steady-state wind field start- 
ing from a set of initial experimental data. Their equations contain 
no time dependent terms. On the other hand, dynamic models 
describe the time evolution of meteorologically variable fields 
starting from an initial state and merging the effects of possible 
time variations of boundary conditions. 

Diagnostic (i.e. kinetic) models follow two main approaches 
[11]: a linearized model which is a simplified steady-state solu- 
tion of the Navier-Stokes equation, and the mass-consistent model 
which analyses available meteorological data, with physical con- 
straints like mass-conservation. It is a general practice to use mass 
consistent models for simulation of stationary three-dimensional 
wind fields in complex terrain. The advantage of a mass consistent 
model compared with primitive equation models is the relatively 
short computing time. However, as a result it is not possible to look 
at complex phenomenon such as turbulence [12]. 

The implementation of a mass consistent wind field is typically 
undertaken using available meteorological data. The procedure has 
been effective for predicting winds for emergency response con- 
taminant dispersion, for weather forecasts, in wind energy studies 
and for storm flooding prediction [13]. Furthermore, some applica- 
tions require that the derived wind field be non divergent. 

The mass-consistent models reconstruct 3-D wind fields 
starting from vertical profiles and near-ground wind measure- 
ments through a two-step procedure [14]. The observed data 
needed for the mass consistent model are provided by an 
interpolation-extrapolation scheme using information available at 
a given site to determine velocity components at each grid point 
above the topography. Typically the available simultaneous data 
consists of several horizontal wind speed and direction measure- 
ments near the earth’s surface, a vertical profile of horizontal wind 
speed and direction and synoptic analysis. 

The wind data interpolation method plays a crucial role in deter- 
mining the final wind field features. Different wind fields can be 
generated starting from the same data and changing only the inter- 
polation method. The adjustment procedure is generally obtained 


Table 1 
Characteristics of the measuring station [19]. 
Site Longitude (°) Latitude (°) Altitude (M) Period of Height of 
measure mast (m) 
(year) 
Ténès 1.33 36.5 17 5 10 


by employing the variational analysis technique [15]. Early research 
was conducted by Sasaki [16-18] who developed the first theo- 
retical model for the calculation of wind speed. Many different 
mass-consistent models have been developed starting from the 
late seventies [19-29]. Various other works relating the results of 
air pollution dispersion modeling in complex terrain are entirely 
dependent on adequate descriptions of the three-dimensional wind 
field [30]. Furthermore, as part of modernization efforts of the 
Atmospheric Release Advisory Capability (ARAC) project, Chan and 
Sugiyama [31] developed a new diagnostic model for generating 
mass consistent wind fields over a continuous terrain. In addition, 
Leone et al. [32] have used wind fields to drive ARAC’s new disper- 
sion model. This model is going to replace the current operational 
code described by Sherman [19] which employs stair-step topogra- 
phy and constant grid spacing. The topography strongly influences 
the boundary layer flow in such a manner that the geographic vari- 
ations either can modify on in some cases generate circulations in 
conjunction with diurnal heating cycles [33,34]. 

The theoretical basis for the mass consistent model was devel- 
oped by Sasaki [16-18]. The general variational analysis formalism 
defines an integral function whose external solution minimizes the 
variance of the difference between the observed and analyzed vari- 
able values subject to physical constraints which are the continuity 
equation and the measured velocities. 

The aim of this paper is to develop a mathematical simulation 
model that can determine the three-dimensional wind velocity 
field over a complex terrain using the Valley of Cheliff in the Tenes 
area in Algeria as a case study. This region has a complex relief with 
a set of mountains and valleys with strong slopes. In the adjustment 
model for wind flow, the relief is analyzed explicitly with the use of 
velocity measurements in the field using a large number of points 
in a three-dimensional grid. The data are interpolated and extrap- 
olated by a least squares method on the grid of an Euler model. The 
numerical results are compared with the measurement values at 
selected points on the real terrain. 


2. Methodology 
2.1. Case study area 


The study area is located in the town of Tenes in Algeria. It is a 
coastal city which is partially surrounded by the Dahra mountain 
range. The topography of this region is complex with a varia- 
tion in altitude between 0 and 700 m (Fig. 1a). The correspondent 
mathematical topographic model is given in Fig. 1b. The data file 
that defines the terrain relief by the three coordinates contains 
374,037 points. Digital processing was used to order the relief struc- 
ture. The site is characterized by an important wind potential area 
[35]. The weather station at this location has 5 years of continu- 
ous measurements of three-hourly speed and wind direction [36]. 
Characteristics of the measuring station are given in Table 1. 


2.2. Mathematical model formulation 


The study is presented as an optimization problem with con- 
straint conditions. The objective function is given by an integral 
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Fig. 1. Digital and mathematical relief of Tenes. 


functional [16-18]: 
E= f(E, v.w, ayavay de (1) 
v 


where F is a function expressed by 


F(u, v, w, A) = a? (u ue)” | a?(v vy | a3 (w wey 


a(g ov w) (2) 
ox Oy Oz 
where u, v, and w are the argument functions which represent 
respectively the components wind velocity in x, y, and z direction of 
the modeling domain V. The mass consistent model constructs an 
observed initial guess of the wind velocities component u®, v? and 
w? over the entire model. The initial estimates are formed from 
an interpolation-extrapolation scheme that uses existing surface 
wind observations contained within the modeled area. The values u, 
v and ware obtained by adjusted scheme from the observed veloc- 
ities by minimizing the functional E. The parameter A(x, y, z) is the 
Lagrange multiplier; and values of œ; are Gauss precision moduli 
taken to be a? =1/ 207. o; is an observation error and/or variances 
of the observed field from the desired adjusted field. Identical Gauss 
precision moduli are assumed for the horizontal directions. The 
ratio w@=(a4/a2)? =(02/04 ) is related to stability and turbulence of 
atmosphere [19]. 

The variational calculus gives the Euler-Lagrange equation for 
u, v and w, respectively [37]: 


oF d (oF), 
ðu dx \ dux ) 


OF d [ OF 

ai 5 
OF d OF 0 

dw dz \ dwy ] 
where 

tie du 

x" OX 

ov 

wy = 5 (4) 
PAN dw 

x Oz 
Eq. (3) are subjected to the boundary conditions: 
W?AVU=0 (5) 


where ñ is the outward positive unit normal to the boundary of vol- 
ume V. Eq. (5) implies that on the boundary, either of the following 
is satisfied: 


VU=0 or A=0 (6) 


Eq. (5) represents appropriate conditions for “flow-through” lateral 
and upper boundaries. Either the normal velocity component vari- 
ation or A must be zero on a boundary. The appropriate condition 
for “flow-through” boundary (upper and lateral boundary) is A =0. 
In this case the derivate of à in the direction across the boundary is 
in general not zero. The lower boundary coincides with the ground 
and represents a “no-flow-through” solid boundary. The observed 
velocity at that boundary is zero [19,25,42]. 
Substituting for F by putting Eq. (2) into Eq. (3), we obtain: 


o, 1a 
u= u" + — 
2q} dx 
p= L9 Fi 
Jaz By n 
pow ue 
= 2a2 dz 


pb E E G) (8) 


The equation for A is derived by differentiating Eq. (7) and sub- 
stituting the results into Eq. (8) giving: 


2 2 2\ 92 0 0 0 
Oh OA as Orr 202 ou po ow (9) 
Ox2  ðy2 a2 ðz2 Ox oy az 

Eq. (9) is an elliptical differential equation (Poisson’s equation). 
It will be solved for A with the specific boundary conditions of Eq. 


(5). The adjusted wind velocity components are calculated using 
Eq. (7). 


2.3. Method of resolution 


This case study can be described as a constraint optimized prob- 
lem where Eq. (1) is the objective function and the constraint 
conditions are the continuity equation and the measured values 
in the volume V. Mass-consistent optimization is used to solve the 
problem. The basic procedure to get the wind field velocities by the 
mass-consistent method is given by the following steps: 

The Lagrange parameter A;;, is determined by Eq. (9) by the 
finite difference method in the discretised volume V which is a rect- 
angular box set on the topographic relief. The grid model is built 
with a choice of Ax and Ay such that the grid takes in account 
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the more interesting topographic detail features of the site (e.g. 
mountains, canyons and valley areas). The computational volume 
is shown in Fig. 1. 

By using a central difference scheme for second derivatives for 
a grid point (i,j,k), Eq. (9) is rewritten as: 


Aiti j,k — 2Àij,k + Ai-1 j,k 
Ax2 


(3 ) Àijik+1 — 2Ai jk + Ài, j,k-1 


Àij+1,k — 2Ai j,k + Aij—1,k 
Ay2 


2 
} a2 Age = —2Q1 £0 (10) 
With the divergence £o of the interpolated wind field is given by: 


ðu) V ðw? 
oy az 


Using central difference for the first derivate £o is written as: 


(11a) 


0 
P Yi gk Yitik 
0 2Ax 


0 0 0 0 
Viet Vijak Wijk Wijke 
2Ay i 2Az 


(11b) 


where the observed values up, es v, pand wp. p Were determined by 


least square weighted method from the measured wind speed. The 
method uses the three nearest measured points to determine the 
observed wind speed at each grid point (ij,k). Once the observed 
velocities are determined, £ọ is introduced in Eq. (10). A set of dif- 
ference equation in Aj; is now obtained at each interior grid point. 
The entire system is solved simultaneously for the value of A; ję. The 
equation set is solved iteratively, using a successive over-relaxation 
method. 

In altitude z the observed values (either measured or adjusted 
values) are extrapolated with the following equation [38]: 


= = z\P 
Ü = Üo ( Z) (12) 
Zo 

where the exponent p is determined by atmospheric stability con- 

ditions or from a least squares fit to multiple level tower data. 
Once the Lagrange parameters are determined, the adjusted 

wind field values are calculated using central difference for Aij 

and mean value for observed wind value of Eq. (7): 


LP 0 o 1 Aistje — Ait jk 
Ui jn gink 2U; j,k uiaj) t Ja? JAx 
1 
1) 0 0 0 1 [Aaja — Ai jk 
Yi jk aij 2v; je tii aad + Fa Dhy (13) 
1 
lo 0 0 1 Ai jk — Aig kt 
Wijk = g Wij + 2Wi jk + Wij) t I JAZ 
2 


The algorithm for calculating the adjusted wind speed values is 
determined as follows: 


e Obtain three-dimensional topographic data of the complex ter- 
rain. 

e Develop three-dimensional discrete reliefs for finite difference 
method. 

èe Determine 3-D numerical box for finite difference method appli- 
cation. 

e Localize measured values of wind velocities on the discrete bot- 
tom boundary of the box. 

e Determine the observed wind velocities (u 
squares method. 

e Calculate the Lagrange parameters in the volume box by finite 
difference method. 

e Use SOR method to resolve set of equations compute 
adjusted velocity components (u, v and w) within the box by 
Euler-Lagrange method (solution of the problem by variational 
form). 


0 v? and w?) by least- 


z Open boundary 


peers, 


Measured value 
(constraint) 


Topographic boundary 


Fig. 2. Finite difference grid of the constraint problem. 
3. Results and discussion 
3.1. Numerical development of the terrain and grid generation 


By letting x, y represent the horizontal direction and z the verti- 
cal direction over a modeling domain V, the data are arranged in the 
y direction and subdivided in columns. The choice of the grid size 
is often determined by the necessity to show the real topographic 
features. Fig. 2 shows a cross section of the computational volume. 
The topography may protrude through the top of the grid poten- 
tially producing unconnected regions within the grid volume. This 
feature can be used to study wind fields in valleys, canyons and 
mountain areas. The numerical relief shown in Fig. 3 indicates that 
the area of interest for this model is a rectangular box set on the 
earth’s surface with the bottom of the box located at the lowest 
topographic point. The dimensions of the box are determined by 
the application requirements and computer storage limitations. 

The computational domain is meshed by 18 x 6 x 10 nodes with 
step Ax = Ay =2000 m and Az=100m, as shown in Fig. 3. The val- 
ues of o; and o were taken to be 1 and 0.01 ms~!, respectively. 
The wind speed and direction at the upper boundary of the grid 
were temporally interpolated from the two vertical profile obser- 
vations and assumed to be spatially constant. The horizontal wind 
measurements were adjusted to a height of 50m above the ter- 
rain using a power law exponent of 0.2 [39]. The constant height 
interpolated and extrapolated wind field is processed to define hor- 
izontal wind components at each grid point about the topography 
boundary. The foregoing applications differ in terrain, meteorolog- 
ical data availability, meteorological condition, grid resolution, and 
convergence criterion. 


3.2. Observed win speed and Lagrange parameters 


The measured values are extracted from results of three- 
dimensional observed velocities (UP ae VP ys WP ad calculated by 
the AILOS software [40]. Fig. 4 shows 11 measured points (black 
symbols) which were chosen on the site at 10 m from the surface. 

The discretised Eq. (9) is solved to give the Lagrange parameters 
values Aj; at each node of the 3D-grid. The matrix coefficients of 
the equation system are diagonally dominant and asymmetric. The 
iterative successive over relaxation, SOR, method [37] was used to 
solve the correspondent equation system. Eq. (9) can be written as: 


Ati = (1 OAT 
1 3 it-1 , 4it-1 ų it-1 \ it-1 i-i it-1 
7 |Aijx PARI HAGI p HARTI + + lambdaitl + BAH + ail )| 
i C 
(14) 
where 
2 

Ai, j,k = 2(a1 AX) £0 (15) 
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Fig. 3. Three-dimensional grid computational volume. 
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Fig. 4. Observed velocities (colored lines) calculated by AILOS software [40] and the 
measured points (11 black symbols). (For interpretation of the references to color 
in this figure legend, the reader is referred to the web version of the article.) 


2 


p=a( $) (16) 
c=ahisa( 2) | (17) 


where Ay = Ax, the coefficients Aij, B and C have known values. 
The variable specifies the current value of the iteration in the SOR 
method. 

The relaxation factor w was determined experimentally and was 
equal to 1.78. This value is chosen so as to obtain the minimum 
number ofiterations. The iterative solution is considered converged 
when the relative change of every Aj; is less than a prescribed 
value. The iterative process is stopped when the condition on the 
solution accuracy is achieved 


jit _ jit-1 


a < Ex (18) 


where ¢€* is a given accuracy value. 
3.3. Wind speed adjusted values 


The numerical technique was applied on the model as shown 
in Fig. 3. The study terrain was very rough, an elevation 


difference of over 500 m between the highest and the lowest layer 
was observed. The model complex relief is composed of a set of 
mountains or hills with a high slope. The major characteristics of 
the flow over a hill were determined not only by the hill shape but 
also by its size and by atmospheric stability. The horizontal wind 
Measurements shown by the colored contour lines in Fig. 4 were 
taken from the observed velocities calculated by AILOS software in 
the study region of Tenes. The accuracy of the wind field is depen- 
dent on the specification of œ and «2. For a stable atmosphere the 
horizontal adjustment were given by a; =1 and œ =10; then the 
adjustment is predominant in the horizontal component instead of 
vertical component since: 


a (2) k (2) 107 (19) 


The developed computer code, applied on the numerical relief 
in shown in Fig. 3, produces nondivergent three-dimensional wind 
at 1080 grid points with the use of a Pentium 4, CPU 2.66 GHz. 


3.3.1. Wind field velocity in horizontal planes 

Fig. 5a shows the horizontal wind field velocity at plane (x, y) 
above the topographic relief. The wind values were adjusted with 
a step height of 100 m in the z direction using Eq. (12). The expo- 
nent p has a value of 0.2 in a period where the atmosphere was 
neutrally steady. The highest wind velocities were observed in the 
intervals [364,366] [352,356] in the x direction. An interesting area 
is a valley with a higher depth from the reference sea level. Fig. 5b 
gives the correspondent wind field which is represented by iso- 
lines of constant velocity. This gives a more precise view of the 
wind field velocity in the horizontal plane (x, y) than the represen- 
tation of velocity given by Fig. 5a. A maximum value of the velocity 
reached 6 m/s in the valley at x = 366 km. It was observed that wind 
direction reversal occurs when the wind shifts across the valley. 
This outcome is also predicted by Finardi et al. [13] for a relief 
with a single valley. The results obtained in our study for differ- 
ent horizontal planes show that the rates decrease gradually as the 
altitude increases and become regular at an altitude of 1000 m. This 
is consistent with the results obtained by Sherman [19]. 


3.3.2. Wind field velocity in vertical planes 

Fig. 6a-e show the velocity distribution in a vertical plane sec- 
tion of the computational volume and in particular the vertical 
motion resulting from terrain features and the wind shear present 
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Fig. 5. Adjusted horizontal wind field in the (x, y) plan above the topography. 


in the vertical profile. The field velocity decreased when the alti- 
tude increase as predicted by Equation (12); the highest velocities 
are located at the wall (flank) of the mountain and the area of the 
valley located in the intervals [364,366] [352,356] as described in 
section of horizontal velocities. This area indicates the location of 
highest power wind. These results are similar to those predicted by 


the two main diagnostic models; MATHEW, developed by Sherman 
[19] and MINERVE, developed by Geai [24] which appear to have 
been the more accurate and the most extensively used. 

Table 2 shows the comparison of wind speeds (m/s) estimated 
using the current model and that estimated using Wasp for the 
Measuring station area Tenes [41]. The current model gave a value 
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Fig. 6. Adjusted wind field velocities in the vertical (x, z) plane. 


Table 2 
Comparison of wind speeds estimated using the current model with that estimated 
by Wasp [41]. 


Coordinate of measurement point (m) Present study Wasp Error 
(m/s) (m/s) [41] 
(X= 3774356, Y=356412.065,0) 3.47 3.69 —6% 


of 3.47 m/s whereas the Wasp method gave a slightly higher value 
of 3.69 m/s [41]. It is noted that the result of the estimate based on 
the current model is lower by 6% compared with that predicted by 
Wasp. 


4. Concluding remark 


A numerical algorithm for the calculation of flow fields was 
developed. The study showed that basic optimization of constraint 
wind flows using a mass-consistent model is an effective method 
for obtaining a 3-D wind simulation over complex terrain. An appli- 
cation was successfully adopted for a complex terrain in the region 
of Tenes in Algeria. The current model gave wind speed values that 
were 6% lower than those estimated by Wasp software. The method 
allowed us to determine the windiest areas in the case study region. 
The precision in modeling complex terrain and the application 
of the scheme difference method on the topographic boundary 
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suggest a significant influence on the wind speed. The model devel- 
oped can also be applied to dispersion of pollutants, and prediction 
of fire spreading. 
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